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ABSTRACT 


With  the  development  of  ocean  surface  remote  sensing,  air-sea  interaction  theory 
and  the  theory  of  underwater  sound  generation  at  the  ocean  surface,  the  potential 
calming  effect  on  surface  gravity  waves  by  raindrop  induced  mixing  has  become  impor¬ 
tant.  The  rain  induced  mixed  layer  was  studied  with  models  based  on  the  turbulent 
kinetic  energy  budget.  A  bulk  mixed  layer  model  proposed  by  Garwood  was  tuned  with 
laboratory  experimental  data  from  Green  and  Houk  (J.  Fluid  Mech.,  1979).  The  tur¬ 
bulent  kinetic  energy  going  into  subsurface  mixing  was  found  to  be  less  than  10%  of  the 
total  raindrop  kinetic  energy.  The  length  scale  for  mixing  is  proportional  to  both 
raindrop  size  and  rain  intensity.  Furthermore,  there  is  some  indication  of  an  initial  pen¬ 
etration  depth  for  raindrops.  Although  the  available  data  was  inadequate  to  complete 
model  development  and  verification,  a  prediction  for  a  hypothetical  situation  in  the 
North  Pacific  is  proposed.  The  diffusion  processes  are  illustrated  by  solving  for  the  dif¬ 
fusion  and  dissipation  terms  of  the  turbulent  kinetic  energy  equation  with  a  finite  dif¬ 
ference  scheme.  New  experiments  are  suggested  to  allow  future  model  development  and 
testing. 
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h  INTRODUCTION 


The  kinetic  energy  of  raindrops  at  the  sea  surface  can  potentially  play  a  role  in  near 
surface  turbulent  mixing.  Such  mixing  will  have  an  effect  on  satellite  remote  sensing, 
air-sea  interaction,  wave  energy  propagation  and  underv^'ater  ambient  noise  production 
through  the  reduction  of  surface  gravity  waves,  temperature  changes  due  to  dilution,  and 
bubble  creation.  Recognizing  this  importance  of  rain  mixing,  it  may  be  useful  to  create 
a  model  which  can  predict  the  turbulent  upper  layer  induced  by  raindrops.  The  results 
might  indicate  the  surface  gravity  wave  damping  rate  which  might  be  used  to  correct 
satellite  data.  As  an  initial  step  a  model  was  built  to  investigate  turbulent  layer  behav¬ 
ior. 

The  energy  of  raindrops  can  go  into  surface  ripples  (Poon  et  al.  1989),  subsurface 
turbulent  mixing  (Chapman  and  Critchlow  1967;  Siscoe  and  Levin  1971),  and  sound 
(Medwin  et  al.  1990),  although  the  partition  of  energy  to  ripples  and  mixing  is  not  yet 
clear.  In  early  studies  about  water  droplets,  the  energy  available  for  mixing  has  been 
closely  related  to  surface  configuration  and  to  drop  shapes.  When  a  water  drop  hits  the 
down  slope  or  trough  of  an  existing  surface  wave  (absorption  effect),  the  drop  energy 
goes  mainly  into  surface  ripples.  On  the  other  hand,  when  the  point  of  initial  contact 
is  upslope  or  at  the  crest  of  a  wave  (reflection  effect),  there  is  a  high  Rayleigh  jet,  and 
most  of  the  energy  goes  to  subsurface  mixing  (Siscoe  and  Levin  1971).  Keedy  (1967) 
shows  that  the  water  drop  changes  its  form  from  spherical  to  vertically  obliated  and 
prolated  while  it  falls.  If  the  impact  occurs  when  the  drop  is  in  a  vertically  obliated  form, 
the  energy  available  for  subsurface  mixing  is  much  larger  than  for  the  case  of  collision 
while  in  a  prolated  form.  These  earlier  observations  indie,  ^  that  some  energy  goes  into 
surface  waves,  some  of  it  goes  into  subsurface  mixing,  and  the  rest  of  it  goes  to  minor 
effects,  e.g.  sound  generation. 

Assuming  that  some  energy  does  go  to  mixing,  Nystuen  (1990)  proposed  e-folding 
times  for  surface  gravity  wave  attenuation  for  various  pairs  of  rain  induced  turbulent 
layer  depth  and  eddy  viscosities  (Fig.  1).  The  most  realistic  values  of  eddy  viscosity  for 
these  predictions  were  not  knowTi.  These  results  show  that  the  strong  mixing  in  a  thin 
turbulent  layer  will  rapidly  damp  short,  less  than  one  meter  wavelength,  gravity  waves. 
In  high  wind  conditions  (larger  background  eddy  viscosity;  curves  G  and  H)  the  in¬ 
creased  attenuation  is  important  only  for  short  wavelengths.  The  short  gravity  waves 


influencing  SAR,  scatterometer,  altimeter  and  passive  radiometers  will  be  damped  if 
substantial  energy  goes  into  mixing. 

Subsurface  mixing  by  rain  on  a  calm  surface  was  investigated  experimentally  by 
Green  and  Houk  (1979)  using  uniform  water  drops  or  a  combination  of  various  drop 
sizes  to  simulate  real  rain.  The  integrated  effects  of  Rayleigh  jets,  vortex  rings  and  sur¬ 
face  gravity  waves  were  assumed  to  influence  entrainment  and  mixed  layer  depth.  The 
experiments  show  that  the  molecular  effects  are  small  compared  with  the  mechanical 
mixing  (inertial  effects)  due  to  the  rain.  The  resulting  mixed  layer  depth  profiles  sug¬ 
gested  that  the  drop  size  has  a  major  effect  on  mixing  and  that  the  intensity  of  rainfall 
is  also  important.  Most  of  the  experiments  were  done  on  fresh  water  for  five  combina¬ 
tions  of  drop  sizes;  fine  spray  where  all  of  the  drop  had  diameters  less  than  1.5  mm,  rain 
with  uniform  drop  sizes  of  2.2,  3.6  and  5.5  mm  diameter  respectively,  and  also  with  a 
mixture  of  the  above  drop  sizes.  Mixed  layer  depth  after  a  fixed  time  interval  was  pro¬ 
portional  to  rain  intensity  and  total  kinetic  energy.  For  larger  drop  sizes,  the  mixing 
was  more  vigorous.  Entrainment  velocity  was  defined  by  the  rate  that  the  mixed  layer 
deepened  with  time.  This  was  inversely  proportional  to  the  bulk  Richardson  number. 
In  salty  water  trials,  mixing  is  less  active  due  to  the  strong  buoyancy  damping  from  the 
salinity  induced  density  difference  between  the  fresh  water  rain  and  the  salty  water. 

In  this  thesis,  these  phenomena  are  studied  with  models.  Two  kinds  of  models  were 
tried.  The  bulk  mixed  layer  model,  originated  by  Kraus  and  Turner  (1967),  is  spatially 
one  dimensional  (time  and  depth)  and  considers  the  integrated  turbulent  kinetic  energy 
in  the  mixed  layer,  assumed  to  be  fully  turbulent.  After  Kraus  and  Turner  (.1967),  a  large 
number  of  bulk  models  have  been  proposed.  These  models  were  designed  to  simulate 
wind  mixing  rather  than  rain  mixing.  Garv'ood  (1977)  presented  an  advanced  bulk 
mixed  layer  model  having  two  original  points.  One  is  the  fraction  of  wind-generated 
turbulent  kinetic  energy  partitioned  to  potential  energy.  This  is  increased  by  means  of 
mixed  layer  deepening,  which  is  dependent  upon  layer  stability.  The  other  is  that  viscous 
dissipation  is  enhanced  for  increased  values  of  the  reciprocal  of  Reynolds  number.  These 
models  are  based  on  the  turbulent  kinetic  energy  (TKE)  equation  derived  from  the 
Navior-Stokes  equation.  Garwood  (1990)  built  a  rain  mixing  model  by  rewriting  the 
NFS  bulk  mixed  layer  model,  adjusting  to  the  case  of  rain  mixing.  Since  real  world  data 
is  almost  non-existent  for  this  problem,  the  experimental  data  of  Green  and  Houk  was 
used  to  tune  this  model.  There  were  four  tuning  constants  which  resulted  in  an  unde¬ 
termined  problem  which  will  be  discussed  in  more  detail. 
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Figure  1.  The  e-folding  Times  for  Surface  Gravity  Wave  Attenuation  as  a  function 
of  wavelength  for  several  choices  of  background  eddy  viscosity,  rain- 
induced  turbulent  layer  viscosity  and  turbulent  layer  thickness.  Curve 
A  shows  the  attenuation  rate  due  to  molecular  mixing.  Curve  B  shows 
the  decrease  in  decay  time  (increased  attenuation  rates)  when  a  layer  of 
weak  mixing  is  present.  Curves  C  -  F  are  for  light  wind  conditions  where 
the  background  eddy  viscosity  is  assumed  to  be  lO'W/s.  Curve  C  has 
no  turbulent  layer;  Curve  D  has  moderate  mixing:  Curves  E  and  F  have 
strong  rain-induced  mixing.  Curves  G  and  H  show  the  much  smaller 
influence  of  a  rain-induced  turbulent  layer  in  moderate  wind  conditions 
when  the  back  ground  eddy  viscosity  is  larger.  Curve  G  is  for  the  case 
of  no  rain  and  Curve  H  has  a  0.1  m  turbulent  layer. 

Another  type  of  model  is  a  diffusion  model  which  encodes  the  diffusion  and  dissi¬ 
pation  terms  of  the  TKE  equation  using  a  finite  difference  scheme.  This  type  of  model 
tries  to  have  the  turbulence  energy  propagate  itself  to  deepen  the  mixed  layer.  Because 
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of  stability  problems  for  the  non-steady  solution,  we  only  used  it  to  show  the  fine 
structure  within  the  mixed  layer  after  Xirst  using  the  bulk  mixed  layer  model  to  predict 
the  mixed  layer  depth. 

The  main  emphasis  for  this  paper  is  to  estimate  the  best  tuning  constant  values  for 
the  rain  model.  This  process  will  be  described  in  Chapter  3  with  a  discussion  of  the 
significance  of  the  chosen  Values  in  Chapter  4.  One  pair  of  optimum  values  will  be  ap¬ 
plied  to  an  hypothetical  ocean  situation  to  provide  a  prediction  for  mixed  layer  forma¬ 
tion  by  rain. 
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11.  EXPERIMENTAL  EVIDENCE 


A.  EXPERIMENTAL  SETUP  OF  GREEN  AND  HOUR 

The  experimental  setup  of  Green  and  Houk  (1979)  is  represented  in  Figure  2.  Rain 
modules,  water-filled  boxes  with  holes  in  the  bottom,  were  used  to  generate  the  artificial 
rain.  The  rain  intensity  was  varied  by  changing  the  water  head  in  the  module.  Each 
module  was  15  cm  deep  and  64.5  cm  square.  Modules  were  constructed  to  produce  ar¬ 
tificial  rain  with  uniform  drop  sizes  (3  different  drop  diameters)  and  to  approximate  na¬ 
tural  rain  with  intensities  of  about  1.25  cm/h  (2%  of  5.5  mm,  13%  of  3.6  mm,  85%  of 
2.2  mm  drops)  and  2.5  cm/h  (4%  of  5.5  mm,  23%  of  3.6  mm,  73%  of  2.2  mm  drops). 
The  experiment  took  place  in  a  vertical  shaft  with  dimensions  of  3  x  6  x  16  m.  The  rain 
module  was  placed  14  m  above  the  receiving  tank.  The  tank  bottom  and  three  of  the 
walls  were  made  of  in.  marine  plywood.  The  fourth  wall  was  made  of  •—  in.  glass. 
One  inch  of  styrofoam  insulation  was  placed  on  the  inside  walls  and  bottom. 

The  temperature  was  measured  by  thermopiles  (22  vertically  aligned  in  the  tank)  and 
sampled  every  14.4  s  with  the  accuracy  of  ±  0.05®C.  A  Fenwall  thermistor,  mounted  on 
a  vertical  rod  (designed  to  move  up  or  down  with  velocity  of  1.1  cm/s)  was  also  used  to 
get  the  vertical  temperature  profile  with  the  same  accuracy.  The  skin  temperature  (top 
10-^  \00  urn)  was  measured  by  a  radiation  thermometer.  Because  the  temperature  of 
rain  falling  through  a  deep  layer  of  air  will  reach  and  remain  at  the  air  wet-bulb  tem¬ 
perature,  the  rain  module  w'as  kept  to  within  0.25‘’C  of  the  wet-bulb  temperature  by  a 
Haake  FK2  constant  temperature  bath.  The  rain  drops  were  generated  using 
hypodermic  needles  to  generate  known  drop  sizes  with  diameters  of  2.2,  3.6  and  5.5  mm. 
The  drop  sizes  were  calibrated  placing  five  drops  into  a  previously  weighed  beaker  of 
olive  oil.  The  rain  intensity  was  determined  by  measuring  the  surface  height  using  a 
plastic  ruler  and  a  low-pass  capacitance  wave  gauge.  In  the  salty  receiving  tank,  the  rain 
intensity  was  measured  by  ruler  only,  and  the  conductivity  gauge  was  used  on  the  mov¬ 
ing  probe,  instead  of  the  thermistor,  to  measure  salinity.  The  receiving  tank  was  left 
undisturbed  for  20  minutes  prior  to  the  start  of  each  experiment.  The  tank  water  tem¬ 
perature  W'as  adjusted  with  ice  prior  to  beginning  the  experiment.  A  plastic  sheet  w'as 
placed  underneath  the  module  before  each  experiment,  and  was  removed  to  start  the 
experiment. 
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Figure  2.  A  Schematic  Diagram  for  the  Experimental  Setup  of  Green  and  Houk. 


Table  1  presents  the  rain  and  environmental  conditions  for  the  experiments  reported 
by  Green  and  Houk  (1979).  The  temperature  of  rain  is  higher  than  that  of  the  water 
surface  in  most  cases;  i.e.,  special  attention  was  paid  to  the  case  of  warm  rain  falling  on 
colder  fresh  water.  The  impact  of  water  drops  creates  a  turbulent  mixed  layer  which 
gradually  deepens  with  time,  mainly  through  the  entrainment  of  the  fairly  quiescent  fluid 
below  the  mixed  layer.  Internal  waves  exist  at  the  thermocline  (the  lower  boundary  of 
the  mixed  layer)  and  capillary-gravity  waves  were  on  the  surface.  The  surface  was  dis¬ 
figured  by  Rayleigh  jets  and  by  many  bursting  bubbles  or  secondary  splashes  associated 
with  the  rain  drop  impacts.  Vortex  rings  are  also  present,  and  it  is  unclear  how  they, 
or  the  waves,  interact  with  the  mixed  layer  turbulence.  Green  and  Houk  studied  the 
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integrated  effects  of  all  of  these  phenomena,  and  dealt  with  these  effects  in  terms  of 
entrainment  and  mixed  layer  depth  (Green  and  Houk  1979). 


Table  1.  CONDITIONS  FOR  THE  FRESH  WATER  RAIN  EXPERIMENT.  All 

_ temperature  are  in  °C _ 


Drop  siie 

Intensity 

Air 

Air 

wet'bulb 

Initial 

surface 

water 

Bottom 

^ater 

(mm) 

(cm/h) 

temperature 

temperature 

temperature 

temperature 

Cl-6 

0-40 

23-1 

16-7 

9-75 

6-25 

OSO 

23-1 

16-7 

8-75 

6-40 

080 

23-1 

16-7 

9-10 

6-60 

090 

23-1 

16-7 

15-30 

13-90 

0'35 

23-7 

16-5 

15-30 

13-70 

090 

23-7 

16'5 

15-35 

13-90 

0'60 

23-7 

Its 

21-15 

21-90 

0'50 

23-7 

16  5 

19-30 

19-90 

2.0 

060 

22-7 

16-5 

8-95 

6-15 

1-50 

22-7 

16-5 

8-90 

6-50 

I'SO 

22-7 

16-3 

9-00 

6-18 

1-80 

22-7 

16-3 

12-40 

10-60 

1*10 

23-0 

15-5 

11-40 

10-00 

1-45 

24-0 

15-7 

11-85 

10-20 

)-40 

24-0 

15-8 

12-50 

10-80 

MO 

24-0 

15-8 

10-15 

8-95 

050 

23-9 

18-0 

16-40 

15-40 

1-60 

25-0 

17-8 

16-10 

15-30 

o-:o 

2S-0 

17-8 

18-90 

19-10 

1-20 

25-0 

17-8 

18-95 

19-15 

3'6 

0-30 

23-3 

150 

7-60 

6-35 

1-45 

23-6 

17-5 

13-25 

11-25 

2-45 

23-9 

17-5 

12-95 

n-40 

3-45 

23-9 

17-5 

13-35 

11-50 

2-20 

23-3 

16-7 

8-95 

6-25 

1-60 

23-3 

16-7 

8-05 

6-60 

3-'!0 

23-3 

16-7 

9-30 

6-60 

1-85 

23-2 

17-2 

15-80 

14-35 

3‘50 

24-3 

17-8 

19-70 

20  00 

2-90 

24-3 

17-8 

21-50 

21-85 

5'i 

0-70 

23-5 

10-7 

9-40 

6  20 

1-20 

23-8 

17-0 

9  30 

6  20 

2-30 

23-8 

17-0 

9-50 

6  60 

2-20 

23-8 

17-0 

12-70 

11-45 

0-40 

23-3 

160 

12-00 

10-75 

MO 

23-3 

160 

12-85 

10-50 

M5 

23-3 

16-0 

13-00 

10  75 

1-40 

23-3 

16-0 

15-10 

13-55 

0-30 

22-8 

150 

14-40 

12-70 

2-10 

22-8 

15-5 

16-95 

17-60 

090 

22-8 

15-5 

16-90 

17-45 

V»riable 

1-20 

23-9 

17-4 

10-25 

6-85 

1-30 

23-9 

17-4 

13-75 

11-40 

1-40 

23-9 

17-4 

16-45 

15-40 

1-40 

23-9 

17-4 

19-00 

19  50 

2-00 

23-9 

17-5 

21-10 

21-60 

2-10 

23-9 

17-5 

12-50 

7-65 

2-55 

23-9 

17-5 

13-10 

11-65 

2-40 

23-6 

17-5 

15-95 

14-35 

B.  DISCUSSION 

It  is  known  that  some  of  the  raindrop  energy  may  go  into  subsurface  motions  not 
associated  with  surface  waves  (Chapman  and  Critchlow  1967).  Siscoe  and  Levin  (1971) 
investigated  the  splashes  that  occur' in  the  presence  of  surface  waves  experimentally. 

They  found  there  are  two  modes;  one  produces  almost  no  Rayleigh  jets,  but  strong  sur¬ 
face  waves,  and  the  other  produces  unusually  high  Rayleigh  jets,  but  less  strong  surface  i 

waves.  The  first  mode  was  defined  to  be  an  “absorption”  event  which  occurs  when  the 
drop  hits  the  trough,  or  the  down  slope  of  a  surface  wave.  The  latter  is  called  a  “re¬ 
flection"  event,  and  occurs  when  the  drop  hits  the  crest  or  the  upslope  of  a  surface  wave. 

Very  little  kinetic  energy  goes  into  the  jet  in  an  absorption  event,  whereas  considerable 
kinetic  energy  goes  into  the  jet  in  a  reflection  event  (Siscoe  and  Levin  1971). 

A  subsurface  vortex  ring,  produced  by  rain  drops,  derives  its  energy  mainly  from  the 
surface  energy  of  the  drops.  There  are  two  conditions  for  optimum  vortex  ring  forma¬ 
tion:  the  drop  must  be  spherical  or  it  must  be  changing  from  an  vertically  oblated  to 
prolate  spheroid  upon  contact  with  the  receiving  tank  water  (Chapman  and  Critchlow 
1967).  Keedy  (1967)  showed  that  the  penetration  of  the  rings  is  related  to  the  initial 
circulation  in  the  drop  and  the  density  difference  of  the  rain  and  receiving  water.  In 
general,  \\1th  larger  circulation  or  smaller  density  difference,  there  is  greater  penetration. 

The  characteristics  of  the  ring  are  highly  dependent  on  the  shape  of  the  drop  when  it  hits 
the  surface.  In  the  case  of  equal  density  of  the  drop  and  the  receiving  water,  the  ring 
diameter  increases  linearly  with  penetration  to  within  about  a  ring  diameter  of  the  crit¬ 
ical  depth.  The  critical  depth  is  the  depth  where  the  rings  stop  (rain  onto  salty  water), 
or  transforms  into  several  smaller  rings  of  very  low  velocity  (on  fresh  water).  For  the 
salty  water  case,  the  ring  diameter  increases  linearly  with  penetration  for  a  time,  but  then 
assumes  a  constant  size  or  even  decreases  with  continued  penetration.  This  suggests 
that  rainfall  on  salty  water  will  not  generate  as  a  deep  a  turbulent  layer  as  rainfall  on 
fresh  water. 

Figure  3  shows  several  vertical  temperature  profiles  obtained  by  Green  and  Houk. 

The  profiles  qualitatively  show  that  mixed  layer  depth  is  closely  related  to  the  drop  size 
and  rain  intensity.  The  initial  temperature  profiles  show  that  the  water  within  6  cm  of 
the  surface  is  slightly  warmer  than  the  underlying  water.  This  shows  that,  before  the  rain 
Stans,  the  sensible  heat  flux  from  the  air  and  the  net  long-wave  radiation  flux  due  to  the 
warm  shaft  walls  more  than  balances  the  cooling  effect  of  evaporation  at  the  surface. 

This  profile  can  be  expected  in  the  real  world  where  a  calm  water  surface  condition  ex¬ 
ists.  Once  rain  drops  start  falling,  mixing  disrupts  this  initial  condition. 
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Figure  3.  Representative  Temperature  Profiles  of  Green  &  Houk's  Experiment. 

(a),(b),(c),(e)  show  thermocouple  profiles  at  t  =  0,  15,  30,  60  min;  (d) 
shows  thermistor  profiles  at  t  =  5,  15,  30,  60  min.  The  rain  intensity 
(upper  number)  and  drop  size  appear  on  each  figure.  The  heat  stored 
near  the  surface  always  increases  with  time. 

Typical  sizes  of  rain  drops  are  from  0.5  mm  to  5.5  mm  in  diameter  (Marshall  and 
Palmer  1948).  Figure  4  shows  the  relation  between  raindrop  diameter  (D),  the  number 
of  drops  per  unit  volume  (No)  and  the  rain  intensity  (R).  depends  on  the  drop  size 
and  rain  intensity.  The  drop  sizes  and  rain  intensities  used  in  Green  and  Houk's  exper¬ 
iments  are  in  Table  1.  The  most  commonly  discussed  conditions  in  this  thesis  have  in- 


tensities  of  0.6  cm/h  for  2.2 mm  drops,  1.6  and  3.7  cm/h  for  3.6  mm  drops  and  2.3  cm/h 
for  .5.5  mm  drops. 


Figure  4.  Typical  Rain  Drop  Sizes. 


The  maximum  vertical  temperature  gradient  increases  with  both  increasing  mixed 
layer  depth  and  time  due  to  the  heat  added  to  the  upper  layer,  and  the  relatively  small 
effect  of  molecular  heat  diffusion  near  the  surface.  The  mixing  by  the  larger  drops  was 
clearly  more  vigorous  than  that  from  the  smaller  drops. 

The  mixed  layer  depth  (H)  was  defined  two  ways.  One  definition  is  to  use  the  mid¬ 
point  of  the  depth  interval  over  which  the  maximum  temperature  gradient  occurred  (this 
definition  tends  to  remove  the  effect  of  molecular  heat  conduction).  The  other  is  the 
level  above  which  90  %  of  the  heat  transferred  through  the  surface  is  stored.  Table  2 
compares  the  mixed  layer  depth  calculated  for  each  of  the  above  two  definitions.  The 
mixed  layer  depth  varies  little  whether  calculated  using  the  heat  storage  definition  or  the 
midpoint  of  temperature  gradient  maximum  definition.  This  suggests  that  the  molecular 
effects  are  small  compared  with  the  mechanical  mixing  by  the  rain. 
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table  2.  COMPARISON  OF  MIXED  LAYER  DEPTHS  (IN  CM)  Calculated 
from  the  heat  storage  (Dh)  and  from  the  maximum  temperature  gradient 
(Dt). _  - 


Drop  size 
(mm) 

Intensity 

(cm/h) 

Dt  (10  min) 

Dh(10 

min) 

Dt  (40  min) 

Dh  (40 
min) 

2.2 

0.6 

5.7 

6.2 

7.2 

8.4 

1.80 

8.0 

8.5 

11.0 

11.2 

1.80 

8.2 

8.9 

11.0 

10.6 

3.6 

1.60 

10.5 

11.7 

14.8 

16.7 

3.70 

15.3 

15.3 

25.0 

24.7 

0.3 

5.0 

7.7 

7.8 

10.2 

5.5 

0.40 

5.5 

5.7 

8.9 

9.8 

1.2 

13.5 

13.9 

27.0 

27.5 

2.30 

28.0 

25.6 

38.5 

37.4  . 

Variable 

1.30 

10.1 

11.3 

16.2 

16.7 

2.10 

19.7 

18.8 

26.0 

25.5 

Figure  5  shows  the  mixed  layer  depth  after  20  min  of  rainfall  into  both  fresh  and  salt 
water.  This  clearly  shows  the  effects  of  drop  size  and  intensity.  The  effects  of  large 
drops  are  especially  evident  in  experiments  in  w'hich  variable  drop  sizes  were  used.  Even 
if  there  are  only  a  few  large  drops,  the  mixed  layer  depth  is  larger  than  for  higher  inten¬ 
sity  rain  experiments  which  contained  only  smaller  drops.  The  mixed  layer  depth  should 
be  strongly  affected  by  the  mechanical  energy  flux  through  the  water  surface; 
O  =  pit)  =  py.  lx  {terminal velocityy  where  1  is  the  rain  intensity,  p  is  the  density  of 
a  rain  drop,  «.  is  the  arbitrary  velocity  scale  used  in  the  turbulence  layer  following  the 
guidelines  of  Green  and  Houk,  and  O  is  the  kinetic  energy  flux  at  the  surface  (Jlm^s). 

Green  and  Houk  defined  the  descent  of  the  thermocline  due  to  the  upward 
entrainment  of  the  lower  fluid  by  an  entrainment  coefficient  E.  This  coefficient  is  de¬ 
fined  as  the  ratio  of  the  descent  rate  to  a  characteristic  velocity  in  the  turbulent 
(ill 

layer  (u.),  i.e.,  E  =  — r-  The  entrainment  coefficient  (E)  is  a  function  of  the  bulk 

(tt  gllAp 

Richardson  number,  which  is  defined  as  Ri  = - The  entrainment  coefficient  was 

piii 

inversely  proportional  to  the  Richardson  number  (slope  is  -1  for  each  drop  size)  regard¬ 
less  of  the  drop  sizes  present  in  the  rain. 
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Figure  5.  Mixed  Layer  Depth  After  20  Minutes  of  Rain.  I  is  the  rain  intensity  in 
cm  and  H  is  the  mixed  layer  depth.  The  empty  circle  symbol  is  for  rain 
with  2.2  mm  drops,  the  empty  box  is  3.6  mm  drops,  'x'  represents  rain 
with  5.5  mnt  drops  and  is  for  rain  with  variable  drop  sizes. 


As  a  summary,  the  Green  and  Houk  experiments  prove  that  larger  rain  drops 
produce  deeper  nrixed  layers  and  suggest  that  a  deeper  mixed  layer  is  generated  in  a  fresh 
water  body  than  that  of  a  corresponding  salt  water  body.  The  entrainment  coefficient 
is  inversely  proportional  to  the  bulk  Richardson  number. 
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III.  MIXED  LAYER  MODELING 

A.  INTRODUCTION  TO  MODEL  TYPES 
1.  The  NPS  Bulk  Mixed  Layer  Model 

The  bulk  mixed  layer  model  is  an  one  dimensional  model.  Typically,  a  fully 
turbulent  mixed  layer  is  assumed  to  be  bounded  above  by  the  air-sea  interface  and  below 
by  a  dynamically  stable  water  mass.  There  is  a  sharp  density  jump  between  this  upper 
layer  and  the  water  below.  This  density  jump  provides  the  conditions  for  Kelvin- 
Helmholtz  instabilities  at  the  interface.  This  region,  which  is  intermittently  turbulent  in 
comparison  with  the  overlying  mixed  layer,  is  called  the  'entrainment  zone'.  An 
entrainment  velocity  usually  is  defined  as  the  rate  of  increase  of  mixed  layer  depth 
(MLD). 

For  a  rain- induced  mixed  layer,  the  kinetic  energy  for  generating  the  turbulence 
comes  from  the  kinetic  energy  of  the  falling  raindrops.  It  is  assumed  that  there  are 
enough  raindrops  to  maintain  the  turbulence  in  the  upper  layer  and  to  deepen  it.  Even 
if  a  large  drop  penetrates  to  a  certain  depth,  the  mixed  layer  produced  by  that  drop 
cannot  be  maintained  due  to  dissipation  and  buoyancy  damping.  If  raindrops  are  sparse 
enough,  the  turbulence  will  disappear  in  a  few  seconds  and  the  upper  layer  of  w'ater  will 
stay  in  a  laminar  state. 

As  the  kinetic  energy  is  transported  from  the  surface  to  the  entrainment  zone, 
some  energy  is  lost  to  dissipation.  This  loss  was  neglected  in  the  prototype  one¬ 
dimensional  model  of  Kraus  and  Turner.  Subsequent  studies  showed  that  this  dissipation 
term  is  important. 

The  concepts  described  above  were  modeled  using  the  turbulent  kinetic  energy 
(TKE)  equation,  vertically  integrated  across  the  mixed  layer.  The  TKE  equation  for 
typical  wind  mixing  is. 


,  .  on 

w'p' 

dt\2  )  1 

u  w  — 
oz 

Po 

■  ^2  Po  ^ 

(A) 

(B) 

(C)  (D) 

(1) 


w'here,  E  =  ^v'^  u,v  and  w  represent  the  velocity  component  of  the  upper  layer, 

p  is  density  and  g  is  gravity.  The  prime  indicates  the  perturbation  term  while  subscript 
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zero  and  overbar  indicate  mean  quantities.  Term  (A)  represents  shear  production  from 
wirid  forcing,  (B)  is  the  buoyancy  flux,  (C)  is  the  TKE  transport  term  and  (D)  is  the 
dissipation  term.  The  wind  stress  source  term  (term  (A))  does  not  exist  in  rain  mixing 
and  raindrop  kinetic  ener^’  flux,  the  (C)  term  evaluated  at  surface  after  vertical  inte¬ 
gration,  takes  the  role  of  the  source  term.  In  the  NFS  rain  induced  bulk  mixed  layer 
model,  hereafter  referred  to  as  the  NFS  bulk  rain  model,  or  more  briefly,  as  the  rain 
model,  the  mixed  layer  is  assumed  to  be  fully  turbulent  as  in  the  other  bulk  models.  In 
the  following  discussion,  the  vertically  integrated  form  of  each  term  of  equation  (1)  will 
come  up. 


2.  Diffusion  Model 

This  type  of  model  is  radically  diflerent  from  the  NFS  bulk  model  in  the  sense 
that  it  attempts  to  have  the  kinetic  energy  of  the  turbulence  diffuse  itself  downward. 
Unlike  the  NFS  bulk  model,  there  is  no  vertical  integration  over  the  mixed  layer  depth. 
Earlier  attempts  with  this  type  of  model  have  had  difficulty  predicting  the  growth  of 
wind-driven  mixed  layers  (Garwood  1990).  Because  of  the  difficulty  in  using  this  model 
to  predict  the  c’-ovt'th  of  the  mixed  layer,  we  attempted  to  use  it  to  investigate  potential 
steady  state  .  ad'tions.  The  purpose  of  this  model  became  mainly  to  show  the  kinetic 
energy  profile  of  the  NFS  bulk  rain  model  more  realistically. 

The  reduced  turbulent  kinetic  energy  equation  is, 


dpk  dpkujf 

6[  dx}( 


^ik 


dxj^  dxj^  ’  6xj( 


)-pc. 


When  we  assume  that  there  are  no  mean  velocities  and  no  horizontal  gradients,  the 
equation  becomes, 


Spk  dp  ^  d 

01  dz  dz 


(pv,-^)-p£ 


i.e.,  the  buoyancy,  diffusion  and  dissipation  terms  remain.  As  a  first  approach  to  mod¬ 
eling  the  TKE  equation,  the  buoyancy  term  was  neglected, 


di 


dz 


/  dk  s 


(2) 


where  the  viscosity  was  v, «  Djk  ,  dissipation  was  c  ==  and  D  was  the  mixing  length 
scale. 
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For  the  stability  analysis,  the  direct  method,  energy  method  or  the  von  Neuman 
method  are  usually  used.  The  direct  method  investigates  the  ratio  of  and  k”  directly 
and  the  criterion  is  that  the  ratio  is  less  than  1.  The  energy  method  tests  v^'hether  the 
sum  of  squares,  is  bounded  or  not.  It  assumes  cyclic  continuity.  The  von 

Neuman  method  assumes  that  the  solution  has  a  wave  form,  substitutes  this  solution 
into  the  finite  difference  equation  and  examines  the  behavior  of  the  amplitude  of  k”  for 
each  wave  number  m.  This  has  to  be  stable  for  all  m  to  be  genuinely  stable.  For  this 
problem,  the  diffusion  equation  is  non-linear.  For  non-linear  cases,  the  energy  method 
usually  works  best,  but  cyclic  continuity  does  not  hold  for  this  situation.  Thus,  for  sim¬ 
plicity,  the  von  Neuman  method  was  chosen  for  this  problem. 

3.  Combining  the  Models 

The  NPS  bulk  rain  model  can  be  used  to  determine  the  MLD  for  various  rain 
mixing  situations,  but  it  does  not  attempt  to  model  the  structure  within  the  mixed  layer 
(i.e.  it  is  vertically  integrated).  The  diffusion  model  explicitly  models  the  diffusion  of 
TKE  coming  through  the  surface  down  to  the  mixed  layer  depth.  It  works  best  with  a 
steady  state  situation  and  does  not  predict  the  MLD  very  well.  Given  these  two  reasons, 
we  attempted  to  combine  the  two  models  to  produce  a  more  realistic  TKE  profile  within 
the  ed  layer. 

Combining  these  two  types  of  model  do  not  include  unifying  the  codes  them¬ 
selves.  The  bulk  rain  model  was  run  first,  and  for  any  particular  time,  the  values  of  TKE 
and  MLD  can  be  determined.  These  two  pieces  of  information  become  the  input  to  the 
diffusion  model.  The  assumption  of  the  NPS  bulk  model  that  there  is  no  TKE  below  the 
entrainment  zone  is  maintained. 


B.  DIFFUSION  MODEL  EFFORT 
1.  Development 

The  modeling  effort  for  equation  (2)  uses  a  finite  difference  scheme.  There  are 
many  different  time  differencing  schemes.  A  list  of  some  of  the  most  basic  and  fre¬ 
quently  used  schemes  are  in  Table  3. 

As  a  first  approach,  we  considered  only  the  diffusion  and  dissipation  terms  in  a 
steady  state  situation. 

Equation  (2)  can  be  scripted  as  follows. 


dk 

dz 


Dk~'2{4!^f+Dki 

oz 


•'2, 

ok 

n  2 
oz 


3 

A-T 
D  • 
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Table  3.  THE  TIME  DIFFERENCING  SCHEMES. 


Name 

Scheme 

Characteristics 

Euler 

/,(n*i)  _  =y(/:W'rtA/) 

Explicit, Truncation  accu¬ 
racy  of  0{Ai) 

Backward 

Implicit,  0{At) 

Trapezoidal 

Implicit,  0[(A/)^] 

Matsuno 

=  Aw-I-Ar/"),  =  ^w+A//''"*)’ , 

/"*')•  -f  1)A0 

Explicit,  O(A0 

Huen 

=  Aw+A//">, 

Explicit,  O[(A0^] 

Leapfrog 

3  level,  O[(A0a 

Adams- 

Bashforth 

)^(.*i)  =  ^w+Af(X/->-/"-i)) 

3  level,  0[(A/)2] 

For  modeling  this  equation  with  a  finite  difierence  scheme,  stability  needs  to  be 
considered.  The  stability  analysis  was  done  using  the  von  Neuman  method.  Except  for 
Adams-Bashforth,  all  of  the  schemes  in  Table  3  were  stable,  The  computer  was  used  to 
chose  the  best  one  based  on  robustness  with  respect  to  the  input  parameters.  A  trial 
was  made  with  a  small  time  period  and  shallow  mixed  layer  depth  to  test  which  scheme 
could  handle  the  widest  range  of  input  parameters.  The  results  suggested  that  the  Huen 
scheme  is  the  best  for  this  problem,  though  it  is  only  conditionally  stable.  The  Huen 
scheme  (in  time),  and  the  centered  scheme,  in  space,  was  chosen  for  this  problem: 

dz  2A2 

,  dk  .2  _ 

d^k 

5-2  ~  a.2 
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'v  "  2 

Therefore  the  finite  difference  form  for  step  1  of  the  Huen  scheme  is, 

= *;+A. 

+ -fr  ( 

and  step  2  becomes, 


The  grid  system  is  chosen  as  shown  in  Fig.  6  because  this  is  a  one-dimensional 
problem. 


2.  Analjlic  Solution 

The  analytic  solution  for  the  steady  state  of  this  case  can  be  derived.  Equation 
(2)  for  steady  state  is 
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Figure  6.  The  Grid  System  for  the  Diffusion  Model. 


oz  oz  D 


=  0. 


Let, 


then  a  coupled  linear  system  is, 


i-rr,  ok 

x,  =  k2D  — 
'  dz 


X2  =  k2 


X 

D 


3xj 
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Solving  this, 


I1.J.  I1.S. 

X,  *=aeV2  D+be  V2  d 

/T~  j-  nr  -  /i-i. 


Since  X2  is  finite  at  z  =  +  00,  a  =  0.  The  solution  is 


D 


Where  it, -(^i)f-/:(z-0). 

The  analytic  solution  shows  that  the  kinetic  energy  profile  depends  on  the  mix¬ 
ing  length  scale  and  the  depth.  With  increasing  depth,  the  kinetic  energy  is  exponentially 
decaying.  If  the  mixing  length  scale  is  larger,  the  kinetic  energy  decreases  more  slowly 
with  depth.  This  solution  was  used  to  help  verify  that  the  diffusion  model  was  working 
properly. 


3.  Stability  Analysis 

For  the  stability  analysis,  the  von  Neuman  method  was  used.  It  is  assumed  that 
the  solution  will  be  of  the  form 


Replacing  the  average  terms 
becomes. 


(4) 


in  equation  (3)  by  kj  ,  equation  (3) 


M. 

dz 


+ i  k;  i. 


(5) 


Substitute  (4)  into  (5),  then 


8k 

dz 


4 


D 

SAz^ 

D 

Az^ 


This  equation  becomes, 


(6) 
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V) 


^  =  (#{  (€“‘^+€-=“'-2)+ 

Using  the  exponential  relations  with  trigonometry,  equation  (7)  becomes, 

(1-  cos^/*Az)+  (1-  cos  /Mz)+  -Jr  }•  (8) 

The  expression  inside  the  brackets  is  a  constant,  C,  thus 

^~-C{k])\.  ■  (?) 

The  first  step  of  the  Huen  scheme  is 

*]«')•  _  *;+4,/ 

The  second  step  is 

Aj''+>>r=A;+-|*Ar(r+/"+‘^‘) 

« ~  Ata"i{l+(1-A/C(;t")‘2')i)]. 

Let  A  -  AtCk^\,  The  condition  for  stability  is  0  ^  .4  ^  1.  The  A  term  is 

{k^)'2{  (1~  cos^'AAz)+  (1-  cos  ikA2)+  ~ 1.  (10) 

If  each  term  of  equation  (10),  has  equal  weighing,  then  each  term  should  be  smaller  than 
-j.  The  first  term, 

gives  the  condition  for  At, 

c(*;)T 

The  second  term. 
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also- gives  a  condition  for  At, 


,A/Z) 


Similarly  for  the  third  term, 


A/:^ 


1_4sL_ 

3  «  i-  ' 

D{kf)2 


Ai< 


1  D 


(12) 


(13) 


Condition  (12)  is  smaller  than  (II)  or  (13)  because  the  value  of  Az  is  usually 
smaller  than  the  value  of  D,  therefore  At  should  satisfy  condition  (12).  In  this  criterion 
for  stability  of  At,  the  TKE  itself  is  included  due  to  the  non-linearity  of  the  problem. 
The  value  for  TKE  changes  depending  on  time  and  depth. 

4.  Initial  Conditions 

The  diffusion  model  was  initially  tried  for  calm  conditions.  The  initial  kinetic 
energy  everj^vhere  within  the  whole  grid  was  set  to  zero.  A  certain  amount  of  kinetic 
energy  flux  was  given  at  the  very  top  grid  point  and  allowed  to  diffuse  downward.  This 
trial  failed  due  to  the  initial  zero  value  of  TKE,  which  appears  in  the  denominator  of 
equation  (3).  To  avoid  this  problem,  another  trial  was  made  which  substituted  the  term 
which  has  k  in  the  denominator  with  a  separate  variable.  These  trials  looked  successful 
but  the  stability  criterion  was  more  severe  than  the  first  trial  forcing  a  limited  range  for 
D,  the  mixing  length  scale.  Because  of  this  limitation,  this  path  was  abandoned  and  in¬ 
stead  a  non-zero  initial  value  for  the  TKE,  in  the  turbulent  layer  was  used.  Eventu¬ 
ally,  the  value  used  for  k^  came  from  the  bulk  model.  The  results  of  the  diffusion  model 
typically  show  a  profile  that  decays  exponentially  initially  with  depth,  but  soon  reaches 
an  equilibrium  and  then  maintains  a  constant  amount  of  TKE  to  the  bottom  of  the  grid. 
Grid  dependence  is  inappropriate  and  thus  the  MLD  value  from  the  bulk  rain  model 
needs  be  used  as  another  input.  Used  in  this  fashion,  this  model  was  useful  for  predict¬ 
ing  the  fine  structure  within  the  mixed  layer.  Further  refinement  of  this  model  might  in¬ 
crease  its  usefulness  (adding  buoyancy,  etc.). 
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C.  THE  NPS  BULK  MODEL 

1.  Review  of  the  Wind  Mixing  Model 

Garwood  (1977)  presented  a  bulk  mixed  layer  model  for  the  upper  ocean  mixed 
layer  due  to  wind  mixing  and  the  buoyancy  flux  of  heat  and  salinity.  These  forcing  terms 
(wind  ihixing  and  buoyancy  fluxes)  have  been  the  principle  focus  of  upper  ocean  layer 
(homogeneous)  modeling  in  recent  decades.  As  an  initial  attempt  to  model  a  mixed  layer 
forced  by  raindrop  kinetic  energy,  the  NPS  bulk  mixed  layer  model  was  adjusted  term 
by  term.  Thus  each  term  of  the  rain-induced  NPS  bulk  mixed  layer  model  will  be  inter¬ 
preted  based  upon  an  understanding  of  the  wind  driven  mixed  layer  model. 

The  depth  of  the  wind-driven  mixed  layer  depends  on  the  balance  between  wind 
forcing,  buoyancy  damping,  usually  represented  by  surface  headng,  and  the  destabiliza¬ 
tion  of  the  lower  interface.  The  resulting  entrainment  is  due  to  Kelvin-Helmholtz  in¬ 
stabilities.  The  best  conditions  for  mixed  layer  development  and  deepening  are  wind 
mixing  while  surface  cooling  is  occurring.  This  combination  is  the  main  cause  for 
latitudinal  and  seasonal  mixed  layer  depth  variations  on  large  scales  and  for  diurnal 
variations  on  small  scales.  The  oceanic  mixed  layer  is  usually  characterized  by  a  uniform 
temperature  layer,  disregarding  the  effect  of  salinity.  The  atmospheric  forcing  terms 
(surface  heat  flux  and  wind  shear)  and  the  vertical  velocity  term  (rate  of  shoaling  or 
deepening  of  the  mixed  layer  bottom)  should  be  given. 

The  bulk  mixed  layer  model  uses  a  vertical  integration  of  equation  (1)  across  the 
turbulent  mixed  layer.  The  source  term  for  wind  mixing  is  the  shear  stress  from  wind 
(term  (A)  of  equation  (1)),  and  the  intermittent  upward  surface  buoyancy  flux  through 
the  surface  (term  (B)  of  equation  (1)).  It  is  usually  assumed  that  a  shear  zone  exists  at 
the  top  of  the  layer  and  that  the  thickness  of  shear  zone,  h„  is  much  shorter  than  the 
mixed  layer  depth.  There  is  assumed  to  be  no  shear  zone  at  the  base  of  the  mixed  layer. 

In  the  entrainment  zone,  the  flux  Richardson  number  {R^  is  proportional  to  R, 
(the  bulk  Richardson  number).  The  flux  Richardson  number  is  defined  as  follows. 


agTw' 


-r~7  du  ,  -r-7  dv 
u  w  — — Hv  w  -r- 

oz  oz 


Thus  term  (A)  from  equation  (1)  can  be  written  as 


Tw' 

Rp 
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The  buoyancy  term,  derived  from  the  equation  of  state  for  sea  water,  is 
P  =  PoU-<T-To)+P{s-So)1 


(14) 


where  p  -  Po+p',  Poi  T'o,  So  are  the  mean  values  of  density,  temperature  and  salinity,  re¬ 
spectively,  and  a  and  /?  are  the  expansion  coefficients.  The  definition  for  buoyancy  is 


b: 


Po 


■g- 


(15) 


By  combining  (14)  and  (15), 


b^g{aM-pAs). 

The  surface  buoyancy  flux  (Bq)  for  the  vertically  integrated  model  is 

Bq  *=  b'w'Q  >=  lagT^o-figT^ol.  (16) 

Ignoring  the  salinity  effect  because  of  the  assumption  that  the  salinity  of  both  upper  and 
lower  layer  is  same  (in  wind  mixing),  the  surface  heat  flux  Tw\  -  -Q  .  Therefore 

Bo  =  -agQ 

and  the  entrainment  buoyancy  flux  (5*)  becomes 

Bh  = 

where  Tw'.^  is  the  heat  flux  at  the  bottom  of  the  mixed  layer. 

The  entrainment  buoyancy  flux  can  be  combined  with  the  shear  stress  term  into 
a  single  expression 

(17) 

Furthermore,  an  entrainment  velocity,  W,,  can  be  defined 


where  h  is  the  MLD  and  w.*  represents  mean  vertical  motion.  The  heat  flux  at  the 
entrainment  zone  can  be  expressed  as. 


23 


(18) 


-Tw'_^  =  ATW, 


Thus,  (17)  becomes 


The  transport  term,  (C),  represents  the  TKE  transport  from  the  surface  to  the 
entrainment  zone.  Some  portion  of  the  energy  will  dissipate  and  rest  will  work  to  deepen 
the  MLD.  This  term  can  be  expressed  as  the  ratio  of  available  TKE,  <  £> ,  to  the  time 
scale,  T  ,  required  to  transport  the  energy.  The  energy  should  be  transported  a  distance 
of  h  (MLD).  Therefore,  the  time  scale  (  x  )  is  taken  to  be  proportional  to  h  divided  by 
the  rms  vertical  velocity  scale,  <  w'^> : 

—2  -L 

X-  a^h<w  >  2. 


Term  (C)  becomes, 

-  5_-L 

<E><  w'  >  2 

_  . 

The  left  hand  side  of  equation  (1),  the  unsteadiness  term, 

<E>  JV, 
expressed  as  aj - - . 

Neglecting  dissipation,  equation  (1)  can  now  be  written  as  follows  in 

entrainment  zone: 


(19) 

( -/j  -(5),  can  be 


^3 


<E>fye 

h 


-a^agATW^+a^ 


<  w' 


rms 


J_  _ 
>2  <  £  > 

1\ 


Where  a,,  a,  and  are  tuning  constants.  If  they  are  assumed  to  be  order  one  as  an  initial 
guess,  then  the  entrainment  velocity  can  be  expressed  as: 


(^)i  <  £  > 

iaghAf+E}  ■ 


(20) 


Thus  the  model  can  be  used  to  predict  the  deepening  of  the  mixed  layer. 


24 


The  dissipation  term  is  the  most  questionable  quantity.  The  net  viscous  dissi¬ 
pation  term  was  entirely  neglected  (Kraus  and  Turner  1967)  or  assumed  to  be  propor¬ 
tional  to  the  source  term  (Niiler  1975).  But  these  models  have  a  flaw  of  a  continual 
increase  in  the  potential  energy  of  density  stratification  and  a  unlimited  deepening  of  the 
mixed  layer  for  no  surface  heat  flux  cases.  The  limiting  value  of  the  mixed  layer  depth 
appears  to  be  proportional  to  the  Monin-Obukhov  length  scale  (L)  in  the  Niiler  model, 
while  the  real  world  situation  it  is  better  estimated  using  the  Ekman  scale.  As  a  more 
advanced  expression,  Zilitinkevich  et  al.  (1979)  assumed  the  bulk  dissipation  to  be 

fA  vlh 

tdz^ia^+a^lfioD—r-, 

where  ).  is  the  Ekman  scale,  v,  is  friction  velocity,  is  the  buoyancy  parameter,  f  is  the 
Coriolis  parameter  and  is  the  stratification  parameter.  This  expression  appears  to  be 
in  better  agreement  with  experimental  data  when  compared  with  earlier  models 
(Zilitinkevich  et  al.  1979). 

The  viscous  dissipation  of  turbulent  energy  is  governed  by  the  cascade  energy 
transfer  from  larger  to  smaller  eddies,  and  occurs  primarily  at  the  length  scale  of  the 
smallest  eddies.  Gant'ood  (1977)  made  an  estimate  of  dissipation  by  taking  the  rate  at 
which  the  largest  eddies  supply  energy  to  the  smaller  eddies  to  be  proportional  to  the 
reciprocal  of  the  time  scale  of  the  largest  eddies.  For  shallow  mixed  layers,  high 
Reynolds  number  an  integral  model  for  dissipation  in  the  mixed  layer,  inde¬ 

pendent  of  viscosity  and  the  small  scale  is 

fo  _  i. 

tdz  =  m^<E>2 

where  E  is  the  vertical  mean  of  turbulent  energy  and  nii  is  a  constant  of  proportionality, 
roughly  equal  to  1  for  w'ind  mixing.  Thus,  c  becomes 

_  ± 
w,  <  £  >  2 

- h - • 

In  the  case  of  deeper  mixed  layers  {R,  ~  1),  the  time  scale  is  the  inverse  of  the 
Coriolis  parameter,  and  simplest  combination  of  the  two  scales  gives 
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Therefore  e  is 


■  <E>l+f<E>  ~ 

- h - _• 

For  the  bulk  mixed  layer  model,  Garwood  chose  the  shallow  mixed  layer  case  (Caspar 
1988)i 

2.  Adapting  the  Bulk  Model  to  Rain  Induced  Mixing 

The  wind  mixing  model  needs  to  be  adjusted  to  the  case  of  rain  induced  mixing. 
For  rain-induced  mixing  there  is  no  mean  shear,  so  term  A  from  equation  (1)  is  zero. 
The  source  term  becomes  the  portion  of  term  (C)  that  is  evaluated  at  the  surface  after 
vertical  integration.  How  much  of  the  raindrop  kinetic  energy  is  available  to  act  as  a 
source  of  turbulence?  When  the  rain  drop  hits  the  water  surface,  subsurface  vortex 
rings,  Rayleigh  jets  and  surface  gravity  waves  are  created.  Only  subsurface  vortex  rings 
represent  energy  that  goes  to  mixing.  This  fraction  of  energy,  >/ ,  is  an  important  tuning 
constant  for  the  source  term  of  the  rain  induced  mixed  layer  model.  The  rest  of  the  drop 
energy  is  presumably  dissipated  at  the  surface  (although  a  tiny  fraction,  10-‘,  goes  into 
sound  energy).  The  source  term  has  the  form: x  y  where  WR  js  the  raindrop 

terminal  velocity. 

There  have  been  many  investigations  of  raindrop  terminal  velocity.  Beard  and 
Prauppacher  (1969)  presented  a  siimmaiy'  of  results  of  these  investigations,  and  collected 
data  which  agreed  quite  closely  with  those  of  Gunn  and  Kinzer  (1949).  Even  though 
there  was  a  flaw  of  using  drops  which  were  allowed  to  fall  in  an  environment  of  air  of 
only  50%  relative  humidity,  the  study  of  Kinzer  and  Gunn  has  been  regarded  as  the 
most  complete  and  widely  quoted  in  the  literature.  In  this  paper  we  use  these  data  to 
assign  terminal  velocities  to  varying  rain  drop  sizes.  The  drop  sizes  used  by  Green  and 
Houk  were  2.2  mm,  3.6  mm  and  5.5  mm  in  diameter.  The  terminal  velocities  for  these 
drop  sizes  are:  690  cm/s  for  a  drop  diameter  of  2.2  mm,  860  cm/s  for  3.6  mm  drops  and 
915  cm/s  for  5.5  mm  drops  (Gunn  and  Kinzer  1949). 

One  of  the  more  different  aspects  of  the  rain  mixing  model  compared  with  the 
wind  mixing  model  is  the  fact  that  salinity  is  a  major  source  of  buoyancy  damping.  For 
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this  reason,  -the  oceanic  mixed  layer  is  defined  as  a  uniform  density  layer  rather  than  a 
uniform  temperature  layer.  From  Equation  (16),  the  surface  buoyancy  working  rate 
becomes 


Bq  =  -pMTp-Tj+pgs^h 

where  the  surface  heat  flux  is T'vv'o  =  —piT^—T),  the  surface  salinity  flux  is  s'w'o  =  /?s  ,  p 
is  the  precipitation  rate  multiplied  by  density  and  divided  by  2,  i.e.  ■—  pi,  and  I  is  rain 
intensity.  The  entrainment  buoyancy  flux  becomes 

Bh  =  WMTi,-D+IIg{s-Si,m. 


Where,  7,  is  precipitation  temperature,  7*  is  bottom  temperature,  s*  is  bottom  salinity, 
T  and  s  are  the  instantaneous  temperature  and  salinity  of  the  mixed  layer,  and  W,  is  the 
entrainment  velocity. 

Keedy's  experiments  intimate  the  fact  that  the  penetration  of  a  rain  drop  is 
mainly  a  vertical  process.  This  makes  it  possible  to  assume  that  the  total  TKE 
£-  When  E  is  substituted  by  in  equation  (20),  the  entrainment  ve¬ 

locity  for  rain  mixing  is 


_  3 

<E>2 

[agM7+£] 


(21) 


where  is  a  tuning  constant. 

As  in  the  wind  mixing  model,  identifying  the  dissipation  term  is  difficult.  For 
rain  mixing,  the  Reynolds  number  (/?,)  is  the  ratio  of  advection  to  viscosity,  which  can 

4ed 

be  expressed  as  R,  -  — ^ —  where  D  is  the  length  scale  of  stirring  by  individual  rain 
drops.  When  the  Reynolds  number  is  large  (unstable),  the  total  dissipation  is  propor¬ 
tional  to  the  total  TKE,  E,  and  the  mixed  layer  depth,  h,  and  is  inversely  proportional 
to  the  mixing  length  scale;  i.e.,  f  tdz  =  M,  h.  If  R,  is  small  (stable),  the  dissipation 

XT—  p 

term  may  be  dominated  the  viscosity,  so  that  / tdz  =  v(  •—  )Vi  =  v  h.  The  dissi¬ 
pation  term  adapted  for  this  model  is  a  combination  of  these  two  cases;  i.e., 

jcdz  =  M,{-^h+-^h).  (22) 
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M,  is  a  tuning  constant.  For  wind  mixing  it  is  order  1  but  for  rain  mixing  its  magnitude 
is  not  known.  Different  values  will  be  allowed  in  this  study. 

The  Obukhov  length  scale  for  the  ocean  nuxed  layer  depth  can  be  represented 
as  the  ratio  of  source  term  to  the  surface  buoyancy  working  rate.  The  precipitation  rate 
cancels  out  when  these  two  terms  are  set  equal.  The  remaining  terms  can  be  expressed 
for  the  mixed  layer  depth  (h)  or  for  ;  i.e., 

Bq  =  -p{a.g{Tp-T)+PgS)h  = 

Source  =  pt](WRf, 


For  h, 


7](WRf 
agiAlTj+figS  ■ 


For  t], 


agAT+ligS 


Therefore,  when  we  know  the  mixed  layer  depth  roughly,  we  can  estimate  tj. 
This  means  that  for  the  steady  state,  the  mixed  layer  depth  is  not  the  function  of  rain 
rate,  but  only  that  of  the  source  and  surface  buoyancy  term.  And  when  the  surface 
buoyancy  works  as  a  production  term  (such  as  cold  rain  on  warm  water),  the  mixed  layer 
depth  can  grow  continuously. 

3.  Tuning  constants 

A  summary  of  the  tuning  constants  used  in  this  model  are  as  follows: 

M^:  dissipation  rate 

M^:  Entrainment  coefficient 

D  :  Mixing  length  scale 

?/  ;  The  fraction  of  input  energy  which  penetrates  the  surface  boundary  layer  and 
works  as  a  source  term 

One  more  possibility  is  the  initial  mixed  layer  depth  1^.  This  term  could  be  used 
to  represent  the  depth  to  which  the  raindrops  penetrate  during  their  initial  impacts.  The 
bulk  model  does  not  require  that  this  be  non-zero,  however  it  was  observed  that  non- 
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zero  values  produced  model  results  closer  to  the  experimental  mixed  layer  depths.  The 

Values  chosen  for  were  arbitrary  (order  centimeters)  or  predicted  from  the  bulk  model 

itself  by  running  the  bulk  rhodel  in  a  spin-up  mode  (very  small  time  step)  until  the  rate 
dk 

of  change  of  TKE,  (■^)»  became  small.  The  bulk  model  was  then  restarted  with  the 
from  the  spin-up  mode  and  with  a  larger  time  step. 

The  dissipation  term  was  composed  in  a  reasonable  manner,  however  there  ex¬ 
ists  the  possibility  that  it  is  proportional  to  the  expression  we  developed.  This  possibility 
is  allowed  through  the  dissipation  coefficient,  A/„  which  can  be  different  from  one  (the 
value  usually  assumed  as  a  first  guess).  Similarly,  the  entrainment  term  has  a  tuning 
constant  coefficient,  M4,  which  is  not  necessarily  unity.  In  the  wind  mixing  model,  these 
coefficients  are  both  order  one.  In  this  rain  model,  we  allowed  both  M,  and  M4  to 
change  as  we  expect  that  the  mixing  due  to  rain  may  occur  on  a  tiny  scale  compared 
with  wind  mixing.  The  entrainment  could  be  smaller,  or  the  dissipation  could  be  larger 
than  with  wind  mixing. 

The  mixing  length  scale,  which  was  assumed  to  be  proportional  to  the  rain  drop 
size,  effects  the  MLD  but  not  very  strongly.  It  was  used  as  a  fine  tuning  constant.  The 
value  for  ri  is  completely  unknown  because  there  are  no  prior  studies  which  predict  what 
portion  of  total  energy  may  go  into  subsurface  mixing. 

The  biggest  obstacle  for  tuning  this  model  was  the  fact  that  there  is  not  enough 
experimental  data  to  be  used  for  tuning.  Most  of  the  data  published  by  Green  and  Houk 
consists  of  only  one  data  point  after  40  minutes  or  two  data  points  (10  min  and  40  min) 
for  a  few  chosen  cases.  This  required  us  to  assume  values  for  several  constants  so  that 
we  could  tune  one  of  them  and  then  retune  the  other  constants  in  an  interactive  manner. 
A/,  and  M^  should  be  universal  constants  for  the  model  given  a  mixing  length  scale  (D) 
and  Ti  can  be  different  in  different  situations. 

4.  Procedure  for  Verification 

The  NPS  bulk  rain  model  used  the  IBM  package  program  TODE'  which  solves 
differential  equations.  Each  term  of  equation  (1)  was  made  into  a  special  function,  con¬ 
stant  or  derivative  so  that  closure  was  satisfied.  Initially  we  assumed  that  the  dissipation 
coefficient  is  one  (M,=  l)  and  that  all  the  rainfall  kinetic  energy  goes  directly  into  sub¬ 
surface  mixing  (tj  -  1).  The  model  was  tuned  to  the  experimental  data  usmg  A/4  at  first 
and  then  fine  tuned  by  varying  D.  Varj’ing  M^  strongly  affected  the  mixed  layer  depth 
while  vaiy'ing  D  only  produced  minor  changes. 
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The  next  trials  were  made  by  setting  M,  and  equal  to  one,  the  values  they 
normally  assume  in  wind  mixing  models,  and  varying  and  the  initial  mixing  depth. 
These  assumptions  implied  that  the  dissipation  and  entrainment  velocity  formulations 
are  the  same  as  in  the  wind  mixing  model,  but  that  the  fraction  of  raindrop  kinetic  en¬ 
ergy  working  to  create  turbulence  could  be  less  than  one.  Tuning  was  done  by  fixing  tj 
so  that  the  MLD  after  10  minutes  matched  the  experimental  data  and  D  had  some  rea¬ 
sonable  value.  The  initial  mixing  depth  was  then  adjusted  to  match  the  40  minute  data. 

The  last  trials  were  made  by  varying  A/j  and  »/  at  the  same  time  with  fixed 
at  a  value  of  one.  The  value  of  D  w'as  set  proportional  to  drop  size  (10  times  drop  size). 
The  initial  mixing  layer  depth,  h^,  was  predicted  from  the  model  in  a  spin-up  mode.  In 


order  to  have  the  model  estimate  an  imtial  value  for  h^,  the  model  was  run  w'lth  a  very 

8E 

short  time  step  (10‘^  sec)  until  -jj-  w'as  small.  The  depth  of  mixed  layer  at  this  time  was 
taken  to  be  The  model  was  restarted  using  a  larger  time  step  with  h^,  the  initial  TKE 


and  the  time  set  to  the  value  given  by  the  spin-up  model.  By  varying  rj  and  A/j,  optimum 


values  for  rj  and  Af,  pairs  were  identified. 
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IV.  RESULTS  AND  DISCUSSION 


■III,  the  following  discussion,  the  results  of  each  set  of  trials  which  were  introduced  in 
the  previous  chapter  will  be  discussed.  Even  though  my  tuning  effort  is  not  a  complete 
one,  I  tried  to  predict  reasonable  sets  of  tuning  constants  for  salty  water  and  fresh  water. 
Additional  experiments  are  needed  for  more  preeise  tuning  of  this  model. 

The  results  of  the  first  trials  (  Af ,  *  1, »/  =  1)  suggested  that  the  mixing  length  scale 
(D)  is  proportional  to  drop  size.  Figure  7  shows  the  relation  of  D  to  drop  size  when  A/4 
is  0.003  in  salty  water.  The  different  values  of  D  for  the  same  drop  size  are  due  to  dif¬ 
ferences  in  rain  intensity.  This  suggested  a  comparison  of  rain  intensity  versus  mixing 
length  scale.  Figure  8  shows  this  comparison.  For  a  given  drop  size,  the  value  of  D 


Figure  7.  Mixing  Length  Scale  (in  cm)  versus  Drop  Sizes  (mm)  in  Salty  Water. 
A/«  is  0.003,  A/,  and  rj  are  equal  to  one. 
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is  proportional  to  the  rain  intensity.  The  relationship  has  a  different  slope  for  each  drop 
size.  For  larger  drop  sizes,  the  slope  is  much  greater  than  for  the  smaller  drop  sizes. 
Thus,  the  mixing  length  scale  is  proportional  to  both  drop  size  and  rain  intensity. 

In  the  first  trials,  M*  was  set  at  3  x  10-^  so  that  the  mixed  layer  depth  of  the  model 
fit  that  of  the  experiments.  However,  there  is  no  physical  reason  why  this  entrainment 
coefficient  should  be  much  smaller  than  in  the  wind  mixing  case  (M*  =  1).  For  this 
reason,  M,  and  Af,  were  fixed  in  the  second  set  of  trials  to  a  value  of  one  and  then  tj  was 
varied.  In  order  to  roughly  reproduce  the  experimental  data,  the  value  of  rj  was  of  order 
10-’. 


RAIN  INTENSITY  (CM/H) 

Figure  8.  The  Behavior  of  Mixing  Length  Scale  versus  Rainfall  Intensity  for  Each 
Drop  Size. 
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Figure  9#  A  Comparison  of  Model  Prediction  to  Experimental  Data.  The  drop  size 
was  2.2  mm  and  the  rain  intensity  was  0.6  cm/h.  The  experimental  data 
is  after  10  and  40  minutes  rainfall. 
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Figure  10.  The  Mixed  Layer  Depth  ProFile  for  Various  Initial  MLDs.  For  3.6  mm 
drop  size  and  1.6  cm,'h  rain  intensity.  The  initial  mixed  layer  depth  was 
chosen  arbitrarily. 
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Green  and  Houk  published  their  experimental  data  at  10  and  40  minutes  only. 
When  I  adjusted  the  model  to  match  to  the  mixed  layer  depth  data  at  40  minutes,  the 
model  at  10  minutes  had  a  smaller  mixed  layer  than  that  of  the  experiments,  i.e.,  the 
model  MLD  value  increases  less  rapidly  than  in  the  laboratory  experiments.  Figure  9 
shows  a  typical  case  where  the  model  is  fitted  to  the  experimental  data  so  that  the  total 
error  (rms  difference)  is  minimized.  This  means  that  the  mixed  layer  depth  needs  to  ar¬ 
rive  at  an  equilibrium  depth  earlier.  This  can  be  accomplished  by  changing  the  initial 
mixed  layer  depth,  .  The  mixed  layer  depth  for  different  initial  mixed  layer  depths  are 
represented  in  Figure  10.  When  the  initial  mixed  layer  depth  is  deeper,  the  MLD 
reached  an  equilibrium  state  sooner.  This  suggests  the  possibility  of  the  existence  of  an 


Table  4.  MLD  FOR  EACH  DROP  SIZE  when  D  and  are  tuned  with  the  condi¬ 
tion  of  A/]  =  Af,  =  1,  and  tj  =  0.002.  The  DT  and  DH  values  are  the  ex- 
_ perimental  data  after  for  10  min./40  min,  of  rain  onto  fresh  water. _ 


Drop  and 
Intensity 

/lo(cm) 

D  (cm) 

h(model) 

DT  (cm) 

DH  (cm) 

2206 

5.0 

0.7 

6.0/8.2 

5.7/7.2 

6.2/8.4 

2218 

B 

0.5 

8.3/10.8 

8.0/11.0 

8.5/11.2 

3603 

5.0 

0.5 

6. 1/8.4 

5.0/7.8 

7.7/10.2 

3616 

7.0 

0.7 

10.3/15.7 

10.5/14.8 

11.7/16.7 

3637 

8.5 

2.0 

15.3/25.7 

15.3/25.0 

15.3/24.7 

5504 

3.5 

0.4 

5.7/9.3 

5.5/8.9 

5.7/9.8 

5512 

4.0 

3.0 

13.0/24.6 

13.5/27.0 

13.9/27.5 

5523 

20.0 

12.0 

25.4/37.3 

28.0/38.5 

25.6/37.4 

V13 

8.5 

1.3 

11.1/16.5 

10.1/16.2 

11.3/16.7 

V21 

15.0 

8.0 

18.7/27.1 

19.7/26.0 

18.8/25.5 
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initial  mixed  layer  depth  {ho).  This  initial  mixed  layer  depth  can  be  interpreted  as  the  rain 
drop  penetration  depth  caused  by  the  initial  physical  impact.  Following  this  suggestion, 
trials  were  made  to  examine  the  behavior  of  the  initial  mixed  layer  depth  and  D. 


Figure  11.  The  Spin-up  Mode  Criterion  for  Rale  of  Change  of  TKE.  Usually  the 
criterion  was  order  of  10*^ . 

Table  4  shows  the  results  for  the  case  of  Mi  =  Mt  —  I,  and  rj  =  0.002.  The  first  two 
digits  in  the  first  column  indicate  the  drop  sizes  (in  mm)  and  the  second  two  digits  indi¬ 
cate  rain  intensity  in  cm/h  without  decimal  point  between  each  two  number.  For  ex¬ 
ample,  3616  means  3.6  mm  diameter  drops  at  a  rainfall  intensity  of  1.6  cm/h.  This 
notation  will  be  used  through  out  this  paper.  V  indicates  variable  drop  sizes  in  the  rain. 
DT  and  DH  are  the  mixed  layer  depth  data  from  Green  and  Houk  where  DT  is  the 
mixed  layer  depth  defined  as  the  depth  of  maximum  temperature  gradient,  while  DH  is 
the  mixed  layer  depth  defined  as  the  depth  at  which  90%  of  the  heat  transferred  through 
the  surface  is  stored.  The  mixed  layer  depths  from  the  model  and  the  experiments  are 
recorded  after  10  minutes  and  40  minutes  of  rainfall.  The  initial  mixing  depth  looks 
roughly  proportional  to  the  drop  size  and  rain  intensity.  Even  though  this  definition  for 


the  .initial  mbdng  depth  is  quite  attractive,  it  introduces  an  additional  free  tuning  con¬ 
stant,  ho ,  which  requires  additional  assumptions. 


Figure  12.  The  Mixed  Layer  Depth  Profile  for  Chosen  Drop  Sizes.  The  number 
5523  means  5.5  mm  drop  diameter  and  2.3  cm/h  rain  intensity.  The 
continuous  data  is  the  model  data  and  the  point  data  is  the  exper¬ 
imental  data  after  10  and  40  minutes  rainfall. 

As  mentioned  earlier,  the  model  does  not  necessarily  need  ho  to  be  treated  as  a  free 
tuning  constant.  Garwood  suggested  a  spin-up  mode  for  the  rain  model.  In  order  to 
have  the  model  estimate  an  initial  mixed  layer  depth,  ho ,  the  model  was  run  with  a  very 
short  time  step  (10"^  sec)  until  was  small  enough.  The  criterion  for  small  enough 
was  on  the  order  of  10-^  For  example.  Figure  11  shows  the  behavior  of  the  change  of 
TKE  ( )  with  time.  Another  run  was  then  started  with  a  larger  time  step  using  the 
initial  time,  ho,  and  HE,  the  total  turbulent  kinetic  energy,  from  the  spin-up  mode.  Using 
the  initial  depth  value  from  the  spin-up  mode,  MLD  profiles  are  presented  in  Figure  12 
for  different  drop  sizes  and  compared  to  the  experimental  data. 


The  last  tuning  trials  focused  on  the  values  for  rj  and  M,  .  The  dissipation  coeffi¬ 
cient,  can  be  one  as  in  the  wnd  mixing  model,  but  might  have  another  value.  The 
value  of  Mi  was  ftxed  to  one  and  the  value  of  D  was  chosen  to  be  10  times  rain  drop 
diameter.  The  value  of  >/  was  varied  from  10-^  to  10-‘  w’hile  the  value  of  M,  was  varied 
from  1  to  100.  For  each  set  of  tj  and  M„  the  mixed  layer  depth  after  10  minutes  and  40 
minutes  data  was  compared  wth  that  of  the  experiments.  Figure  13  shows  the  optimum 
pair  of  A/,  and  tj  for  each  selected  drop  size.  The  optimum  pair  was  defined  as  the  M, 
and  t]  pair  which  has  the  minimum  least  square  error  when  compared  to  the  exper¬ 
imental  MLD  after  10  and  40  minutes.  As  A/,  is  believed  to  be  a  universal  constant, 
these  results  suggest  that  tj  is  proportional  to  drop  size  and  rain  intensity.  Table  5  shows 
the  combined  (summation  over  all  drop  sizes)  least  square  error  for  each  pair  of  A/,  and 
ri  from  Fig.  3.  A  band  of  optimal  Af,  and  tj  pairs  is  apparent. 

We  don't  expect  the  value  of  Af,  to  be  much  larger  than  that  of  the  wind  mixing  case 
(Af,  =  1).  Therefore,  tj  is  probably  less  than  0.01  implying  that  the  TKE  used  in  subsur¬ 
face  mixmg  is  less  than  1%  of  the  total  kinetic  energy  of  the  rain.  If  the  dissipation  rate 
is  ver\'  large  (Af,  ~  100),  then  rj  may  be  as  high  as  10%. 
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Table  5.  COMBINED  LEAST  SQUARE  ERROR  FOR  EACH  M,  AND  >1  PAIR. 


Af, 

n 

0.001 

0.002 

0.003 

-0.005 

0.008 

0.01 

0.03 

0.05 

0.08 

1 

162.5 

551.5 

1177.8 

2630.9 

2 

163.3 

5 

712.9 

302.4 

174.8 

164.3 

388.5 

604.7 

10 

422.1 

214.6 

141.8 

162.9 

1221.4 

15 

621.7 

366.1 

195.8 

151.0 

611.6 

20 

502.6 

288.1 

211.0 

361.5 

936.0 

■ 

25 

379.1 

286.9 

241.9 

619.6 

30 

360.7 

204.2 

430.1 

40 

178.5 

239 

50 

187.4 

162.2 

60 

241.3 

139.7 

70 

141.6 

201.5 

80 

157.0 

162.8 

90 

182.3 

144.2 

100 

208.4 

138.2 

To  demonstrate  combining  the  diffusion  model  with  the  bulk  model,  a  pair  of  tj  and 
M,  from  Table  5  was  chosen  for  furthf"  analysis.  While  this  pair  (>?  =  0.08,  M,  =  100) 
had  the  minimum  combined  least  square  error  in  Table  5,  these  values  for  rj  and  A/,  may 
not  be  the  true  values  for  these  constants.  Nevertheless,  the  bulk  rain  model  was  run 
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with  these  constants  for  20  minutes  of  rainfall.  After  20  minutes  of  rainfall,  the  mixed 
layer  depth  and  total  TKE  from  the  bulk  model  were  recorded  and  used  as  the  initial 
conditions  for  the  diffusion  model.  Figure  14  presents  the  TKE  profile  from  the  dif¬ 
fusion  model  for  four  different  rainfall  conditions.  The  results  suggest  that  rain  con¬ 
taining  larger  drop  sizes  generates  deeper  mixed  layer  depths  and  more  turbulent  mixing. 


Figure  14.  The  TKE  Profile  from  the  Diffusion  Model.  For  chosen  drop  sizes  in 
fresh  water.  The  depth  is  in  cm.  The  notation  2206  means  that  the 
raindrop  size  is  2.2  mm  and  the  rain  intensity  is  0.6  cm/h.  The  other 
numbers  have  similar  interpretation. 


The  Richardson  number  reflects  stability.  The  entrainment  velocities  were  studied 
with  respect  to  this  number  in  the  experiments  of  Green  and  Houk.  Their  results  show 
that  the  entrainment  velocity  is  proportional  to  the  reciprocal  (slope  is  -I)  of  the 
Richardson  number.  That  the  entrainment  velocity  is  bigger  in  the  unstable  state  than 
in  the  stable  state  appears  to  be  a  common  phenomena.  The  NFS  bulk  rain  model  was 
investigated  for  this  relationship. 
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Figure  15.  The  Richardson  Number  Versus  Entrainment  Velocity 

For  the  rain  model,  the  relation  between  entrainment  velocity  (WE)  and  the  Richardson 
number  is  presented  in  Figure  15.  This  result  shows  that  the  entrainment  velocity  is 
roughly  proportional  to  the  reciprocal  of  the  Richardson  number  for  the  bulk  rain 
model. 

The  most  helpful  function  of  a  model  may  be  prediction.  With  this  in  mind,  I  chose 
a  Ml,  ri  pair  from  the  previous  trial  (M,=  100,  »7=0.08).  We  assume  a  hypothetical  rain 
with  a  uniform  raindrop  size  of  roughly  3.6  mm  and  an  intensity  of  1.6  cm,'h  in  the 
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Figure  16.  The  Time  History  of  Each  Parameter  in  the  Rain  Bulk  Model:  (A)js 

for  temperature  (T)  and  salinity  (S),  (B)  is  for  buoyancy  working  rate 
at  the  surface  (BO)  and  at  the  bottom  of  MLD  (BH),  (C)  is  for  mixed 
layer  depth  (H)  and  total  TKE  (HE)  and  (D)  is  for  the  TKE  (E)  and 
entrainment  velocity(WE). 
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is  chosen  to  be  34.4  ppt,  an  average  value  at  mid-latitudes  in  the  north  Pacific  (Pickard 
and  jEmery  1982).  The  spin  up  time  was  chosen  to  be  15  seconds.  After  this  spin  up 
time,  the  initial  mixed  layer  depth  was  5.4  cm.  Figure  16  shows  the  time  history  for  each 
parameter  in  the  bulk  model.  As  the  figure  shows,  the  salinity  within  the  mixed  layer, 
S,  decreases  continuously  with  time  due  to  dilution.  The  temperature,  T,  quickly  reaches 
an  equilibrium  state  and  maintains  that  level.  The  buoyancy  flux  at  the  surface,  BO, 
decays  wth  time  due  to  the  mixing.  The  buoyancy  flux  at  the  bottom  of  entrainment 
zone,  BH,  increases  rapidly  initially,  but  then  slowly  increases.  The  entrainment  veloc¬ 
ity,  WE,  decreases  with  time  and  the  average  kinetic  energy,  E,  also  decays  with  time. 
However,  the  total  kinetic  energy,  HE,  increases  because  the  mixed  layer  depth  is  in¬ 
creasing  continually. 


Finally,  after  5,  15,  30  and  60  minutes  of  rain,  values  for  MLD  and  TKE  were  re¬ 
corded  and  used  as  inputs  for  the  diffusion  model.  The  mixed  layer  TKE  profiles  from 
the  diffusion  model  are  show-n  in  Fig.  17. 
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V.  CONCLUSIONS  AND  RECOMMENDATIONS 

There  is  clear  evidence  that  rainfall  generates  a  thin  turbulent  mixed  layer  at  the 
ocean  surface.  Such  a  mixed  layer  will  effect  air-sea  interaction  processes,  however  the 
strength  and  character  of  that  turbulent  layer  is  unknown. 

This  paper  has  attenipted  to  model  the  formation  of  the  rain-induced  turbulent  layer 
by  modifying  the  existing  NFS  wind  mixing  bulk  model  to  the  situation  of  rain  mixing. 
The  resulting  NFS  rain  mixing  bulk  model  has  several  tuning  constants  associated  with 
it:  ri,  the  percent  of  the  raindrop  kinetic  energy  available  for  mixing;  A/i,the  dissipation 
rate;  M4 ,  the  entrainment  rate;  D,  the  mixing  scale  and,  perhaps  ,  h^,  the  initial  pene¬ 
tration  depth  of  the  rain  drops.  By  examining  the  values  that  these  tuning  constants 
were  required  to  assume  to  have  the  model  match  experimental  data,  we  can  draw  some 
conclusions  about  the  rain-induced  mixed  layer. 

For  the  NFS  wind  mixing  model,  the  dissipation  rate  and  entrainment  rate  constants 
are  order  one  (i.e.  =  1).  If  the  rain  model  also  assumes  that  they  are  order  one, 

then  j/  is  order  10-^  This  implies  that  only  1%  of  the  kinetic  energy  of  the  raindrops 
will  go  into  subsurface  mixing.  The  rest  of  the  energy  stays  at  the  surface  in  the  form 
of  ripples  or  surface  waves.  Alternatively  M4  may  be  small  (order  10-’),  which  seems 
physically  unsatisfactory,  or  can  be  larger  (order  10  ~  100).  In  fact,  given  a  value  for 
Mj,  7]  has  an  optimum  value  which  may  be  proportional  to  drop  size  and  rain  intensity. 
However,  with  A/,  =  100,  rj  is  still  only  0.08.  Therefore  the  conclusion  that  very'  little 
of  the  raindrop  kinetic  energy  is  available  for  subsurface  mixing  seems  very'  solid. 

Another  conclusion  is  that  the  mixing  length  scale  is  proportional  to  both  drop  size 
and  rainfall  intensity.  There  is  also  some  indication  that  there  is  an  initial  penetration 
depth  which  is  order  of  centimeters  and  also  proportional  to  both  drop  size  and  rainfall 
intensity. 

An  attempt  was  made  to  model  the  rain  induced  mixed  layer  with  a  diffusion  model. 
In  this  study  the  diffusion  model  took  a  trivial  role  of  providing  internal  fine  structure 
given  initial  conditions  from  the  NFS  rain  mixing  bulk  model.  There  is  possibility  of 
further  development  for  this  model;  in  particular,  the  addition  of  the  buoyancy  terms 
should  be  attempted. 

Finally  a  prediction  for  the  formation  of  a  rain-induced  mixed  layer  is  provided  given 
some  hypothetical  rain  conditions  in  the  North  Facific  Ocean. 
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The  experimental  data  used  to  constrain  the  model  consisted  of  two  mixed  layer 
depths  (after  10  and  40  minutes)  for  a  variety  of  artificial  laboratory  rainfall  situations. 
This  w'as  not  sufficient  to  fix  the  values  of  the  tuning  constants  used  in  the  model. 
Further  experimental  data  should  remedy  this  problem.  In  particular  it  should  be  pos¬ 
sible  to  design  laboratory  experimental  procedures  to  identify  tj,  the  energy  available  for 
mixing;  D,  the  mixing  length  scale;  h^,  the  initial  penetration  depth  of  the  raindrops  and 
M«,  the  entrainment  rate.  To  monitor  the  formation  of  the  mixed  layer,  the  data  should 
be  recorded  continuously,  rather  than  at  just  two  discrete  time  points. 
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APPENDIX  A.  THE  NPS  BULK  RAIN  INDUCED  MIXED  LAYER 

MODEL 


VARIABLES  &  INITIAL  CONDITIONS: 

HE  =  .0 

H  =  .  lOOOOOOOOOD-01 
S  =  34.40000000 
T  =  20. 00000000 
TIME  =  . 0 

CONSTANTS: 

AG  =  . 2000000000 
BG  =  . 8000000000 
SB  =  34.40000000 
TP  =  18. 00000000 
TB  =  19. 00000000 
Ml  =  100.0000000 
M4  =  1. 000000000 
D  =  3.600000000 
NU  =  . lOOOOOOOOOD-01 
WR  =  860. 0000000 
P  =  .5000000000D-03 
EBS  =  .8000000000D-01 

SPECIAL  FUNCTIONS: 

E  =  HE/H 

WE  =  M4*E*SQRT(E)/(E+H*(AG*(T-TB)+BG’*f(SB-S))) 
BO  =  -P*(AG*(TP-T)+BG*S)*H 
BH  =  WE*H*(AG*(TB-T)+BG*(S-SB)) 

DIS  =  -2.*M1*(NU*HE/D**2+E*SQRT(E)*H/D) 

DHE  =  SOURCE+BO+BH+DIS 
SOURCE  =  EBS*P*WR**2 

DERIVATIVES: 

D(HE  /D(TIME  )  =  = 

SOURCE+BH+BO+DIS 
D(H  /D(TIME  )  =  = 

WE+P 

D(S  /D(TIME  )  =  = 

(WE*(SB-S)-P*S)/H 
D(T  /D(TIME  )  =  = 

(WE*(TB-T)+P*(TP-T))/H 

OUTPUTS: 

TITLE:  RAIN  MIXING 


TABULATE:  TIME 

H 

HE 

E  WE 

DHE 

SOURCE 

AT  INTERVAL 

. lOOOOOOOOOD-01 

PLOT:  H  HE 

DHE 

DIS 

AGAINST:  TIME 

AT  INTERVAL 

1. 000000000 

END  CALCULATION  WHEN  TIME  .  GE.  30. 0000 


DIS 
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APPENDIX  B.  CODE  OF  THE  DIFFUSION  MODEL 


*  * 

*  TURBULENT  KINETIC  ENERGY  DIFFUSION  MODEL  * 

*  * 

*<c^***i»r*<r<r*)fe*'A’<i*A'*'***'A'**A'<f**<iic*A'**A'<r*'<f/r*"A"/r*'**<r<r*'A"A'**********<r^ 
A  * 


*  THIS  MODEL  DEALT  THE  DIFFUSION  OF  TKE  WITH  THE  GIVEN  * 

*  INITIAL  TKE  VALUE  AND  MIXED  LAYER  DEPTH.  THE  UNIT  OF  MI-  * 

*  XED  LAYER  DEPTH  IS  CM,  AND  OF  RAIN  INTENSITY  IS  CM/H.  * 

A  A 


A  AAAAA'AAAAAAAAAA'AAAAAAAA'AAAAAA'A'AAAAAAAAAAAAAAAAAAAAAAAAAAArA  A 
A  A 


A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 


TKE  :  VARIABLE  WHICH  IS  REPRESENTING  THE  TURBULENT  KIN-  * 

ETIG  ENERGY  (TKE).  * 

TKEO  :  VARIABLE  WHICH  IS  REPRESENTING  THE  TKE  VALUE  OF  * 

ONE  TIME  STEP  BEFORE.  * 

TKEN  ;  VARIABLE  WHICH  IS  REPRESENTING  THE  TKE  VALUE  OF  * 

ONE  TIME  STEP  AFTER.  * 

TKEO  ;  THE  VALUE  OF  TKE  EXIST  INITIALLY.  * 

TKET  ;  TiiE  VALUE  OF  TOTAL  TKE  AT  A  GIVEN  TIME.  GIVEN  BY  * 

USER.  * 

DT  J  TIME  STEP.  * 

DZ  :  SPACE  STEP.  * 

RHO  J  DENSITY  OF  WATER  IN  CGS  UNIT  (*  1.0).  * 

WR  :  THE  TERMINAL  VELOCITY  OF  RAIN  DROP.  * 

RR  :  THE  VALUE  RAIN  RATE.  THE  UNIT  IS  CM/H.  * 

D  :  MIXING  LENGTH  SCALE  IN  CM.  * 

ETA  ;  THE  PORTION  OF  INPUT  ENERGY  WHICH  IS  USED  FOR  * 

SUBSURFACE  MIXING  AS  A  SOURCE  TERM.  * 

DEP  ;  DENSITY  OF  WATER  IN  CGS  UNIT  (=  1.0).  ‘  * 

P  :  SPACE  STEP  WHICH  USER  WANT  TO  PRINT  OUT.  * 

* 


c 

c  Define  the  maximum  time  and  space  value, 
c 

PARAMETER( JM=1000 , JMP 1=1000 ,NPTS=1000 ) 
c 

c  Define  variables, 
c 

REAL*8  TKEO(0; JMP1),TKEN(0;  JMP1),TKE(JM) 

REAL*8  DT, DZ, WR, RHO, RR,D, ETA, DEP, P, TKEO, TKET 
c 

COMMON  JMAXP 
c 

c  Prepare  the  output  data  file, 
c 

CALL  EXCMSCFILEDEF  1  DISK  THESIS  DATA  A’) 
c 

c  Read  in  data  of  mixed  layer  depth,  total  tke,  mixing  length  scale. 
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rain  rate,  and  raindrop  terminal  velocity. 

PRINT*,  ’WHAT  IS  THE  VALUE  OF  MIXED  LAYER  DEPTH?' 

READ(5,*)  JMAX 

PRINT*,  ’what  is  the  VALUE  OF  TURBULENT  KINETIC  ENERGY  AT  THE 
1  SURFACE? ' 

READ(5,*)  TKET 

PRINT*,  'what  is  the  VALUE  OF  MIXING  LENGTH  SCALE?’ 

READ(5,*)  D 

PRINT*,  'what  is  the  VALUE  OF  RAIN  DROP  TERMINAL  VELOCITY?' 
READ(5,*)  WR 

PRINT*,  ’what  is  the  RAIN  INTENSITY?’ 

READ(5,*)  RR 

Set  the  values  of  constants. 

JMAXP  =  JMAX  +  1 

DT  =  0. lD-02 

DZ  =  0. IDOO 

RH0=  l.ODOO 

ETA=  l.OD-03 

P  =  RR*RHO/(7.2D+03) 

TKEO  =  P*WR**2 


Set  the  initial  condition. 

DO  11  J=1,JMAX 
TKE(J)  =0.0 
CONTINUE 


Set  the  value  at  time=0 


DO  22  J=1,JMAX 
TKEO(J)  =  TKET 
22  CONTINUE 


c  Set  the  boundary  condition  using  the  subroutine  'BNDRY'. 
c 

CALL  BNDRY('.TKEO,ETA) 
c 

c  Main  loop  (time  loop) 
c 

DO  100  N=1,NPTS 

c  set  the  boundary  conditions. 

CALL  BNDRY(TKEN,ETA) 

c  compute  Tke  using  huen  scheme. 

CALL  HUEN( TKEO, TKEN,DT,DZ,D,NPTS, JMAX, JMAXP) 
c  tansfer  the  tke  for  next  time  loop. 

DO  55  1=1, JMAX 

TKEO(I)  =  TKEN(I) 

TKE(I)  =  TKEO(I) 

55  CONTINUE 

c  print  out  the  results  to  data  file  using 

c  format. 

IF  (M0D(N,1000).EQ.0)  THEN 
C  WRITEd,*)  'SEC=',N/1000 
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DO  111  J=l,:tMAX/3 
DEP  =  -3*J*DZ 
WRITE(1,222)  DEP,TKEC3*J) 
777  CONTINUE 

END  IF 

222  FORMATC  ',6D12.3) 
iOO  CONTINUE 
999  STOP 
END 


Subroutine  BNDRY 


SUBROUTINE  BNDRY(KO,ETA) 


REALMS  KO(0:JMAXP) 
REAL*8  ETA 

COMMON  JMAXP 

K0(0)  =  ETA*K0 
KO(JMAX+l)  =0.0 

RETURN 

END 


*  * 

*  Subroutine  HUEN  * 

*  * 
irMritIfitirkiritIfificifidfkiticicirkiticirititirk'Hifkirk'ltirirkit'kitit 


SUBROUTINE  HUEN(KO,KN,DT,DZ,D,NPTS, JM) 
c 

REAL*8  K0(0:  JMAXP) ,KN(0: JMAXP) ,DT,DZ,D,ETA 
INTEGER  JMAXP 
c 

COMMON  JMAXP 
C 

c  Huen  scheme 

DO  33  J=1,JM 

KN(J)  =  K0(J)+D*DT*(K0(J+1)-K0(J-1))**2/ 

1  ( 8 . *SQRT( ( KO( J+1 ) +KO( J- 1 ) ) / 2 .  DOO ) *DZ**2 ) 

2  +  D*DT*SQRT((K0(J+1)+K0(J-1))/2.D00)* 

3  (KO(J+l)-2.*KO(J)+KO(J-l))/DZ**2  - 

4  DT*SQRT(((K0(J+1)+K0(J-1))/2.D00)**3)/d 
33  CONTINUE 


DO  99  J=1,JM 

KN(J)  =  KO(J)  +  (DT/(16.*DZ**2))*(D*((KN(J+1). 
5  KN(J-1))**2/SQRT((KN(J+1)+KN(J-1))/2.D00) 
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+  ( K0( J+1 ) -K0( J- 1 ) )**2/SQRT( C  K0( J+1 )+K0( J- 1) ) / 2 . DOO ) ) 

+  8.D00*D*(SQRT((KN(J+1)+KN(J-1))/2.D00)*(KN(J+1) 

-  2.*KN(J)+KN(J-1))+SQRT((KO(J+1)+KO(J-1))/2.DOO) 
*(K0(J+1)-2.D00*K0(J)+K0(J-1)))  -  8.*DZ**2* 
(SQRT((CKN(J+l)+KN(J-l))/2.  D00)**3)+K0( J)*SQRT(((KO( J+1) 
+KO(J-1))/2.DOO)**3))) 

WRITE(6,*)  'J=',J,'K(J)=',KN(J) 

CONTINUE 


RETURN 

EITO 


APPENDIX  C.  FLOW  CHART  FOR  THE  DIFFUSION  MODEL 


Main  Program 
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APPENDIX  D.  FLOW  CHART  FOR  THE  DIFFUSION  MODEL 

(SUBROUTINE) 


Subroutine  HUEN 
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